Quantum algorithm for exact Monte Carlo sampling 
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We build a quantum algorithm which uses the Grover quantum search procedure in order to 
sample the exact equilibrium distribution of a wide range of classical statistical mechanics systems. 
The algorithm is based on recently developed exact Monte Carlo sampling methods, and yields a 
polynomial gain compared to classical procedures. 
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The possibility of using quantum mechanics to treat in- 
formation and perform computation has attracted great 
interest in the recent past (see e.g. [Ij] for a review). 
Quantum algorithms have been devised, which solve com- 
putational problems faster than their classical counter- 
part, such as the factorization algorithm of Shor 
However, relatively few problems have been identified 
which are amenable to quantum speedup. While many 
works have proposed methods to simulate quantum sys- 
tems using a quantum processor (see e.g. [l| and ref- 
erences therein), fewer have tried to build quantum al- 
gorithms to speed up classical physical problems {3j. In 
particular, statistical physics is the source of many com- 
putational problems which have led to great efforts to 
develop efficient classical algorithms. For example, the 
goal of many Monte Carlo algorithms is to sample a con- 
figuration set D, from an equilibrium probability distribu- 
tion TT [4]. It is therefore important to explore the pos- 
sibilities to use quantum computers in order to speed up 
such problems. Some quantum algorithms have been pro- 
posed to approximate the partition functions of certain 
statistical physics models or even obtain it exactly 
in very specific cases @. In the very recent past, many 
works have focused on quantum algorithms implementing 
classical Markov Chain Monte Carlo (MCMC) methods 
through quantum walks (7l-[To|. In general, these methods 
give a quadratic gain compared to classical simulations. 

Here we consider another type of MCMC algorithm 
recently developed, the "coupling from the past" (CFTP) 
procedure of Propp and Wilson, which leads in finite time 
to the exact equilibrium distribution [11]. We propose a 
quantum algorithm combining this CFTP procedure and 
the quantum search procedure of Grover [12| , enabling a 
quadratic speedup over the classical algorithm, without 
using quantum walks. Our method enables to sample 
the exact equilibrium distribution in finite time for a wide 
class of systems, while previous algorithms either provide 
an approximate version whose error has to be controlled 
[5, 7-10], or are restricted to specific models [^j. Our 
algorithm is also rather simple compared to these other 
methods, while yielding a comparable polynomial gain. 

One of the key issues in classical MCMC algorithms 
is that they must be iterated sufficiently many times so 
that the final state is a "typical" configuration, in other 



words has a probability distribution close to the station- 
ary one, TT, independently of the initial state. In order 
to get close to the correct distribution tt, one should be 
able to know when sufficient convergence is achieved. In 
a few particular cases, it is possible to calculate analyt- 
ically the relaxation time of the algorithm, i.e. the typ- 
ical time needed to reach stationarity. But in practice, 
estimating or bounding relaxation times is a notoriously 
difficult mathematical problem (T3l - [l8| and one has to 
rely on heuristic arguments to infer that stationarity has 
approximately been reached. 

An elegant alternative way to circumvent this issue has 
been proposed by Propp and Wilson in 1996 flT', 19']. As 
detailed below, the CFTP technique is a reformulation 
of the MCMC procedure that generates exact samples, 
in the sense that they are exactly distributed according 
to the stationary distribution tt. Thus successive calls 
of the algorithm generate totally uncorrelated samples 
(see below). The basic idea is to run the Markov chain 
"from the past" , from a time —t up to t = 0. Now 
suppose that there exists a time —T such that at t = all 
states have coalesced (or "coupled"), i.e. their evolution 
through the algorithm has led to the same state Xc of the 
configuration set fi. Then any initial configuration at t — 
—00 would lead to the same state Xc, which can thus be 
seen as the result of an infinite time simulation. The state 
Xc is consequently distributed exactly according to the 
stationary distribution. The difficulty of the procedure 
dwells in the necessity to track the evolution of the whole 
set fl, whereas in standard Monte Carlo sampling only 
one state of H. is tracked. 

When stored in a computer memory, the configuration 
set is always finite. Thus we will consider a discrete time 
Markov chain [2^ on a finite configuration set £7 of car- 
dinality N. Our quantum algorithm consists in replacing 
the classical evolution of the N states of 57 by a quantum 
evolution of a single quantum state, namely the superpo- 
sition of the N ones, \x); then the Grover quan- 
tum search procedure is applied on top of this quantum 
evolution to find efficiently if the system has coalesced. 

Let us now detail the classical CFTP algorithm. With- 
out loss of generality, a state x G fl (with f2 of cardinal N) 
can be coded by n classical bits 6^ = 0, 1 as a; 61 . . . 6„ 
and iV < 2". Let P be the transition matrix of the 
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Markov chain. Its elements are the transition proba- 
bihties P{x,y), with P{x,y) > being the conditional 
probability that the chain is in the state x at time t + 1 
given that it was in the state y at time t. The chain 
is supposed to be reversible, which means that it satis- 
fies the detailed balance condition 0, HO] : there exists a 
probability distribution on fl, denoted by tt, such that 
TT{x)P{y,x) = ■n{y)P{x,y) for all states x and y. This 
condition ensures that tt is a stationary distribution. We 
assume that tt exists and is unique, in which case it coin- 
cides with the equilibrium distribution (see [l^l for fur- 
ther details). If P{x,T\xq,Q) is the probability that the 
chain is in the state x at time T given that it was in the 
state xq aX t = 0, then [20] 

lim F(a;,r|a;o,0) = tt{x). (1) 

T^oo 

This result is central in traditional Monte Carlo sam- 
pling: if the algorithm (the Markov chain) is iterated 
long enough, then the probability distribution of its final 
state is close to the stationary one. 

A Monte Carlo step at time t can be seen as a map 
/t : 51 — ^ il, determined by a randomly generated pa- 
rameter at as /t(-) = (/)(•,«(). Thus once the random 
numbers at are set, each step of the algorithm is com- 
pletely deterministic. If T is the duration (number of 
steps), then the algorithm is entirely coded by the map 
Ft = /t o • ■ • o /2 ° /i ■ A standard way of perform- 
ing Monte Carlo sampling consists in following the dy- 
namics of a single initial state during a sufficiently large 
time and averaging a physical observable O over time 
iterates. In this case, the statistical error on the nu- 
merical measure of (O) can be estimated using the re- 
laxation time To of O. Indeed, the algorithm behaves 
as if roughly T /tq independent realizations were mea- 
sured, leading to an error Erro — \J2to/TI\0, where 
T is the total simulation duration and AO the standard 
deviation of O [3]. In principle, tq can itself be numer- 
ically measured through the correlation function of O, 
Co{s) = {0{t + s)0{t))t - {Of oc exp(-s/Tci). How- 
ever, Co{s) is itself an equilibrium quantity that can be 
measured only if the system has reached stationarity, and 
the measured to may be only representative of a long 
transient regime instead of equilibrium. This is partic- 
ularly critical in disordered, glassy systems where it is 
impossible to ascertain whether equilibrium has indeed 
been reached, because the system is likely to get trapped 
for a long time in the many metastable states |2l|-|23j . 

Instead of iterating the states in the future as in the 
traditional method explained just above, the CFTP pro- 
cedure goes from the past: we suppose that we have at 
our disposal a sequence of random numbers q;_i, a_2, . . . 
before starting the algorithm. The CFTP algorithm con- 
structs the iterates of all the states a; G H until they have 
all reached the same state ("coalescence"). The essential 
subroutine of the algorithm [let us call it n(T)] calculates 
the N computational paths from time —T up to time 
through the map Gt = f-io. . - of-T and tests at the end 



whether all histories have coalesced to the same state. If 
the coalescence test fails, the same procedure is started 
again from an earlier time. The algorithm reads: 

T = 
repeat 

T = T + AT (go AT steps back in time) 
n(T) (follow all N paths and test coalescence) 
until coalescence is achieved 

At coalescence, T is such that one has Gt{x) = Xc 
for any x & fl. Thus for any t' < —T , G-ti{x) = Xc- 
In particular, limf'_)._oo G-f (a;) = Xc- Therefore Xc can 
be seen as the result of a Monte Carlo algorithm of in- 
finite duration and is exactly distributed according to 
TT. It is proven that with probability 1 the algorithm re- 
turns a value in finite time [ll|. Successive calls of the 
algorithm return totally uncorrelated samples, and (O) 
is now calculated by averaging over realizations instead 
of time. The statistical error is now perfectly controlled: 
Erro = where R is the number of realizations (inde- 
pendent calls of the algorithm). 

The CFTP can be applied to any MCMC problem (H]. 
In specific instances, it is sufficient to follow the history 
of a small subset of the N states. This is the case e.g. if 
a partial order on fl makes it sufficient to follow only 
extremal configurations. Unfortunately, in general the 
N states should be followed in parallel, which represents 
an often prohibitive computational cost. We will propose 
below a quantum algorithm reducing this cost. 

We denote the average coalescence time of the algo- 
rithm by T. To what extent is it related to relaxation 
times To as discussed above? In principle, to depends 
on the observable O. But To is bounded above by (and 
in general on the same order of magnitude as) the re- 
laxation time of the Markov chain, denoted by r [25j . 
This latter time measures the speed of convergence of 
the probability distribution to tt and is equal to the in- 
verse of the first gap of the transition matrix P [3 ■ 
Furthermore, r is itself bounded above by (and in general 
on the same order of magnitude as) f [17], which makes 
CFTP-type techniques so useful to estimate convergence 
rates, even on theoretical grounds. All in all, generally 
speaking, f ~ r ^ tq. Running the algorithm yields an 
accurate estimate of f and thus of relaxation times. 

Let us now turn to our quantum algorithm. The es- 
sential subroutine n(r) of the classical CFTP algorithm 
follows the history of the N states x G fl and tests 
coalescence of all states. The quantum algorithm will 
start from a register holding the N initial states \x), 
< X < N — 1, coded on n qubits, and follow each 
history in parallel. A second register holds the results of 
the successive application of the maps ft . Calculation of 
the iterates is done with the help of ancilla registers. To 
illustrate the computational scheme, we will first special- 
ize our presentation to the case of the Ising model with 
Glauber dynamics. In this case the observable O could 
be e.g. the magnetization. Starting from the totally sep- 
arable state, Hadamard gates are applied to each qubit to 
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put the register into an equal superposition of all states. 
Suppose that after t iterations the quantum state reads 



\xm,Tixm\om, 



(2) 



x=Q 



where Ht T = /-T+t-i ° • • f-T- The next Monte Carlo 
step consists in flipping a certain spin i with some proba- 
bility function of the energy difference AE between con- 
figuration Ht^T{x) and the same configuration with spin 
i fiipped. This spin i may be chosen at random; alter- 
natively one can consider each spin one after another 
(sequential sweep). The spin is flipped with probabil- 
ity p = [1 -I- exp(/3Ai?)]~-'^ or left unchanged with prob- 
ability 1 — p. Here /3 is the inverse temperature. In 
terms of quantum registers, one has to evaluate the en- 
ergies associated with configurations |y) — \Ht^T{x)) and 
= X^'^ly) (where the Pauh matrix X^*^ flips spin i), 
by arithmetic operations controlled by the second regis- 
ter. The probability p is then calculated on the third reg- 
ister and the random number at, uniformly distributed 
on [0,1], is put on the fourth register. The sign of at — p 
is computed, and the one-qubit fifth register \s) is set to 
|0) if a — p > and |1) if a — p < 0. The bit-fiip matrix 
X^*) is applied, controlled by \s). The last three registers 
are then reset to |0) in the usual way by running the op- 
erations backwards. The circuit in Fig. [T] shows one step 
of the iteration algorithm before reset of these registers. 



|0) 
|0) 
|0) 



j4„ 



x) 
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sign(p - a)) 



FIG. 1: Circuit for one step of the Monte-Carlo algorithm. 
The unitary operator Up calculates probability p, Va is the 
operator implementing the random numbers a, and Ap-a is 
a modified adder circuit that gives the sign of p — a. 



After T steps the quantum state reads 



E n)|Gt(x)). 



(3) 



If T is such that all Gt{x) are equal then the histories 
have coalesced and a measure of the second register yields 
an element in fl distributed exactly according to the sta- 
tionary distribution tt. Since all iterations are performed 
in parallel this step requires an average time f. The cru- 
cial point in the CFTP algorithm is that the exact dis- 
tribution is obtained if and only if all Gt{x) are equal. 
Consider for instance the case where the Gt{x) take two 
different values, say j/i and ?/2- Then one might detect 
that the states have not coalesced onto a unique value 
as soon as one obtains different results after measuring 



the second register upon repeated runs of the procedure 
(with the same random numbers) . If the states have not 
coalesced then the process has to be restarted from an 
earlier time. Obviously, only in the case where the proba- 
bilities of measuring yi and ?/2 sue both high will different 
outcomes be obtained quickly upon measurement. In the 
extreme case where there is a unique xq € such that 
Gt{xo) = 1/2 while all other x verify Gt{x) = yi, almost 
all measurements of the second register will give yi , and 
the state will be almost indistinguishable from the state 
where all Gt{x) are equal. Since the CFTP algorithm 
requires to distinguish these cases, the idea is to use the 
Grover algorithm to amplify the probability amplitude of 
the unknown noncoalesced state |a;o), so that it can be de- 
tected. If we first measure the second register and consis- 
tently get the value yi then our aim is to detect whether 
all x verify Gt{x) = yi or not. Since we know the value of 
yi , we can attach to our quantum state a one-qubit regis- 
ter in the state \z) — -^(|0) — |1)). We then perform the 
operation \x)\Gt{x))\z) i-7> \x)\Gt{x))\z + ip{x)) , where 
ip{x) = if Gt{x) — j/i and '^(x) = 1 otherwise, and 
erase all the registers but the first one. This gives a phase 
e*'^ to states which do not verify Gt{x) = yi- One can 
thus apply Grover iterations to magnify the amplitude 
of the noncoalesced states. The whole sequence above 
corresponds to one "oracle" step of the Grover iteration, 
and has to be performed using the same random numbers 
at- 

What is the speed-up on n(r) obtained by this proce- 
dure? Suppose that after T time steps M states x are not 
coalesced. We consider the case M <^ N, since otherwise 
noncoalescence is easily detected by a few measurements 
or even classically. Then 0{y^ N/M) Grover iterations 
are required (even though M is unknown [lEl), each us- 
ing 0(T) operations. The total number of quantum gates 
in this case is Ty/N/M. As the random numbers can 
be generated classically once and for all beforehand, their 
computational cost is ~ T and thus negligible. 

As explained above, the complete algorithm proceeds 
by performing a certain number of calls of n(r) until co- 
alescence. Let us evaluate the speedup of the quantum 
algorithm for n(r) compared to the classical procedure. 
We consider that preliminary runs enable quickly to es- 
timate a suitable AT < f which is gradually improved. 
To iterate classically the Monte Carlo steps on one com- 
putational path will cost a certain number of computa- 
tional operations g{N). Let us first suppose that the 
dynamics is rapid, with a short coalescence time f with 
g{N) ^ \n'^ N. This means that the dynamics is poly- 
nomial in the physical system size n. The classical algo- 
rithm will need to follow ~ N/M computational paths 
to detect the ones which did not coalesce. The total cost 
is g{N)N/M ~ (TV In" N)/M. In contrast, in the case 
of the quantum algorithm, 0{^J N /M) calls to the oracle 
are required, and the oracle performs the TV evolutions 
in parallel in also ~ g{N) ~ ln° N operations. The total 
speedup of the quantum algorithm will be 0{y/ N/M), 
a quadratic gain. If the dynamics is torpid, with a long 
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coalescence time f with g{N) ^ N'^, then the classical 
algorithm requires ^ N'^^^/M operations whereas the 
quantum one costs only ^ N'^^^^^ /^/M elementary steps. 
The relative gain gets smaller for increasing c, going from 
quadratic for small c to almost zero for c — > cx). 

The estimates above assume that all initial states, 
stored in a n-bit register and coded by all numbers x 
between and N — 1, are physically admissible. This is 
the case e.g. for spin problems. However, in many inter- 
esting cases one must inspect a large number of states 
{N = 2"), of which only a small subset {Na = 
6 < 1) are admissible. This is the case e.g. for dimer, 
spanning tree or hard core lattice gas problems. In this 
situation, one can use a modified version of the above 
quantum algorithm. Indeed, starting from an equal su- 
perposition of the N states, one can build in the Grover 
oracle [before each n(r) step] a subroutine which recog- 
nizes the nonadmissible states, and overwrites in this case 
the second register with a known admissible state Xa , so 
that the latter is used as an initial state in the dynam- 
ics. In this case, the quantum algorithm for M = 0(1) 
still requires ~ g{N)^/N steps, since the Grover search 
is applied on the whole Hilbert space of dimension N. 
To obtain both the equilibrium distribution and the re- 
laxation time, the classical algorithm needs N oper- 
ations to identify admissible states, and ^ g{N)N'' op- 
erations to run the CFTP on the admissible states. If 



g{N) ~ ln° A'', the gain is unchanged. If g{N) ~ iV°, it is 
max(c-|-6, l)/(c-|-l/2). The gain is polynomial in all cases 
except for 6 < 1/2 and c > 1/2 (very long relaxation time 
and very few admissible states). 

The algorithm proposed here presents several advan- 
tages compared to the recently proposed method for sim- 
ulating Markov Chain systems 7HlOl|. These procedures 
use quantum walks to approximate the stationary distri- 
bution, in a time typically quadratically faster than the 
classical convergence time. The speedup over classical 
computation is therefore comparable, but, in our case, 
we obtain a sampling of the exact stationary distribution 
rather than an approximate version of it with errors that 
have to be controlled. Our quantum algorithm is also 
simpler. Another advantage of our method is that the 
relaxation time is directly related to the average coales- 
cence time [l3] and thus can be accurately measured. 

It has been proven that the Grover algorithm is opti- 
mal, in the sense that the number of calls to the oracle 
cannot be lower (see [ij and references therein). Our 
quantum algorithm can therefore be improved only by 
speeding up the oracle part. This may be possible by 
combining our algorithm with techniques used in the 
other approaches of (5l-[l0l|. As the algorithm developed 
here is very general, applying to a wide class of systems 
without any structure taken into account, tailored algo- 
rithms may achieve a larger gain in specific cases. 
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